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ABSTRACT 


B method 1s proposeGst Of eppemescimation Of the transite: 
function of a linear, time-invariant system with no numera- 
tor dynamics, from random ineuu and output data. The mevmod 
employed utilizes the Fast Fourier Transform Algorithm and 
Least Squares Estimation to obtain the coefficients of the 
system's transfer function. 

The procedure has been modeled in FORTRAN IV on an 
IBM-360 computer. The results of simulation show the 
feasibility of estimating the order of the transfer function 


and its €@errrerears. 
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IT. INTRODWUSTION 


The problem of identifying a system (or plant) has been 
the object of much thougnvt and. eCirort ase rociCce oy ene 
Tdentif®eation in Automatic Control Systems symoosium (|29 
The effort has been primarily expended on methods which are 
Suitable to a time domain analysis. Identification in the 
frequency domain has in general involved lengthy computa- 
tions to obtain a Fourier Transform if the computations are 
done digitally, or a tradeoff between many parallel narrow 
band filters or repetitive playback of the system input 
and output if the analysis is done by analog methods. The 
method to be used in this investigation has been constrained 
to make use of the Fast Fourier Transform algorithm. There- 
fore digital computations and techniques logically follow. 
The Fast Fourier Transform (FFT) algorithm reduces the 
number of complex multiplications required in the evaluation 
of the Fourier Transform and thus may provide a means of 
performing system identification on-line or approximately in 


real time. The FFT is discussed in Appendix A. 


[ins SBEROBLEM STATEMENT 


The mathematical model of a linear, time-invariant 
system may be expressed as a differential equation or as 
a transfer function. The determination of the model 
parameters irom obsenvatiens son the systcmewielh dopiiowbe 
modeled is generally called an identification problem. 
ThesGussinne st maith on weeriionadeies sieleamode li niga © feomien Simone 
which is characterized bya vector differential equation, 
such that the coefficients of the equivalent transfer 
function Ca@iabe estimated. 


Let a linear system transfer function be written as 


oy = prey > Cf) = ud d (jw) . (1) 
Tne. equimzalkeont qn order difflerenfial equation is w rigeen 
; a ata) = v(t) . (2) 
n=0 at 
dy = 1.0 


From (2) the equivalent N-vector differential equation 


(3) Can de oObtamned,. 





GER eo tT kes O- 
to i 
‘Ko(t) 3 10 0 il a 0 X(t) re 
= ee) vaaie) 
| Pi) 
ce gO CO 1 | ey y(t) | 0 
Le : 
Ry) | dy -9y -dy Teen) syle) | Lee 
(3) 
OR 


mG) = Ax(t) + By(e) 
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MUev=-m loll bimakeher ven? (3) ae 


G 
x(t) = @(t,to)x(tg) +f (t,t) + Br v(r)ar, (4) 


0 
where x(t) is the system state vector and the element x, (t) 
is the system output x(t). The discrete version of (4), 
Which is derived in Appendix B, is used as the equivalent 
ce aa) o 


are estimated. The estimation scheme is represented by 


representation of (1) and the coefficients ad 


Figure l. 
The investigation is to be constrained so that v(t) is 
broad-band white noise and the Fast Fourier Transform 


algorithm must be used in the identification scheme. 


alike 


chiles Bares Ie 
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Estimation Scheme Representation 


Tiles. SYSTEM, MODEL BGR. sie As COME DAEs 


A linear, time-invariant system represented by the 
vector daitrerential equation (5)"hasmeemermioa Gas its 
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(6) 


If v(t) is deterministic and the convolution integral 
can be solved, then a discrete solution of (5) can be 
A(t-ty) 
modeled on a digital computer. e¢ can be converted 
to a square matrix by using the Laplace transform method 
[2, pg. 316]. However, since v(t) is not deterministic the 
convolution integral cannot be solved, so an approximation 


moO (6) is required. 


Considering only the unforced response of (6), 


A(t-ty) 
x(t) =e x(ty) . (7) 
At some instant of time t = a = vO Teel , 
_ _AMT 


which is identical to the unforced response of equation (3) 


of Appendix B where x(MT) = x(t, + MT) = x(ty), 


Ss 


CMe) =se. xm) (9) 
Cee (0) =e non 

m(oT n= ex T) =e anc) 

x(MT) =e" x(t)) (10) 


Therefore, since (10) is equivalent to (8) which is the 


discrete solution @Gieg@/) > then (9) is the’ discrete solution 


for the unforced response of (6) for t = ha ty + (mt1)T 
rr in the m@bation of Meoemciease. bo=— (mie) Ts 
The convolution integral from (6), 
ie = 
g(t) = f htt area en C Tia) 
0 


cannot be forced to yield an exact discrete solution without 
performing the indicated integration. Therefore, the sam- 
pled data approach in Appendix B is necessary. Normally 
this approach yields 


2 cA(T-Temtyar, Cis 


estate) 12) ec) 
0 
where v(mT) is assumed to be constant for mT<t<(m+1)T 
[2, pg. 340). This approach presupposes that the bandwidth 
of v(t) is small compared to = Since v(t) is assumed to 
be broadband noise, T should be small. However, as T 


decreases, the number of iterations required in the simula- 


tion for a given time span of x(t) increases. Each iteration 


14 


requires that both the unforced and forced responses of 
equation (6) be computed. T should then be large to min- 
imize computation and small to allow the maximum transfer 
of information to the system concerning the variations of 
v(t). Equation (12) of Appendix B was chosen as a compro- 
mise. By using this equation, v(t) can be sampled at a 

il 


rate greater than 7 and the unforced response of equation 


(6) need only be computed once every T time interval. 


ike, 


LV. SYSTEM TDENTIFPICATION SCHEME 


A continuous, time=-invamilteant system is to, besadentigied 
by estimating the coefficients,of ats transfer function. 
Thesinputatoathews ys GeiaiSepeesc— band whitesmolcemof Kniewn 
mean and variance. The identification scheme is depicted 


by Figure 2, wheresktadenotes, the Foumier Transformai.e .4 


i Jil ic) 

X(f me F [oat 

v(t) - broadband white noise 
x(t) - system output 


Mies osjccmmoneulls Investigation is To make useuem 
the Fast Fourier Transform (FFT, Appendix A) in place of F, 
and leastesGUuameoeBscimationv in the compulavions 1or tie 


transier humerlom cOoctTi:LCLlents . 
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In Figure 2, vit) 18 the oulgplP of 2 icandom process 
In particular it is the output of 7a Baussian random precess. 
Since v(t) is random, then x(t), V(f) and X(f) are random 
processes. V(f) and X(f) are not themselves deterministic; 


however, from Goldman, [3, pe. 279] 
H(t) = XG) gigas, (13) 


where H(f) is the frequency response of the unknown system. 


mimmce H(t) is essumed toewee Of the form 


= eee lime 
H(f) = + eae 
Sd (j2nf) 


then (13) can be written as [Appendix C, equation (2)] 
Gig) es DIC) Baia (14) 
From (14) a residue r(f) is formed, 
BD wa eae DP je—ey (Tt), ny) 


Through the Fourier Transform, measurements on x(t) and 
v(t) yield X(f) and V(f). <A modified least squares esti- 
mation procedure then uses r(f) from (15) and its complex 


conjugate to form 


F 
i aaa ap aaa ewig) Oi, (kG) 


ef 


where r*(f) is the complex conjugate of r(f). The object 
is then to minimize Q@ with respect to the qd, > the 
Ggietant cents of D(f). 

Since the FFT [Appendix A] is used in place of the 
Fourier Transform of Figure 2 and the unknown system is 
modeled in accordance with Appendix B, the identification 
scheme takes on the form of Figure 3. 

Since u(t) is broadband noise, the low-pass filter is 
a necessary addition to the identification scheme to prevent 
foldover of the alias frequencies of V;- The low-pass 
filter used is a discrete version of an R-C filter with 
its break point on a Bode plot being greater than the 
expected bandwidth of the system under investigation. The 
sampling rate of the inputs to the FRT's was then chosen 
to be at least twice as great as the low-pass filter cut- 
off frequency, that frequency at which the filter output 
drops by 3db. 

X. and Wy are the results of the Discrete Fourier Trans- 
formew@mecay) “and uy) Wiichwere observed through a time 


window of t seconds, i.e. 


Ps 
NW 


a 
X(f,) = XC) 


= 
Hl 

It 
= 
lie 


Acre, 


where t is the time interval over which u(t) and x(t) are 
observed. As the index i goes from -tW to tW, X. and V; 
represent the complex values of the Fourier Transforms of 


x(t) and u(t) at the discrete frequencies -W to +W. 
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The X,'s and an exist only at frequencies which are 
multiples of =. If x(t) and v(t) contain only frequencies 
which are exact multiples of = then the X,'s and V,'Ss 
will be the exact Fourier coefficients of x(t) and v(t) [6]. 
However, since it is unlikely that the only frequencies 
present are multiples of =, these other frequency components 
will @ause the X.'s and V,'s to be altered. It would be 
advisable to weight the frequency spectra of x(t) and 
v(t) wath the@ilanmine for hanning funetions [4, 5, 6, Appen— 
dix A] to minimize the effects of these other frequency 
components. This weighting was not performed, although it 
is admittedly advisable, since computer computation time 
would be increased. The trade-off of less computation time 
for increased accuracy was felt to be justified since much 
of the computer time used in this investigation had to be 


cnareeer tre ovner projects. 
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V, RESULTS 


Many of the variables in the computer simulation of the 
"unknown" plant and the identification scheme [Table Dl, 
Appendix D] are peculiar to the approach taken in this 
investigation. For an actual application of the identifi- 
cation scheme the "unknown" plant is not simulated and only 
the sampling rate and number of samples to be taken must 
be determined. The expected or known signal bandwidth 
dictates the sampling rate, where the minimum rate is the 
Nyquist sampling rate. An increase in the number of 
samples taken increases the frequency resolution of the 
FFT output and the accuracy of the estimated transfer 
function coefficients. However, it also increases compu- 
tation time and computer memory required. 

A tenth order system such as equation (17) was assumed 
rorgeaeh trial; 

Ges) at 0) 
Gy 1 oe 


Be od} 
=0 


lie 


al 


The low-pass filter shown in Figure 3 is represented by 


eeviavion (ld). 


<< 


CS eet CO 
U(S)  S+27- 100 (18) 


Table 1 presents the most encouraging of the results 


Sf rims Mavestigavion. Trials 1] and 2 represeny firsv=o7eon 


eal 


systems with real S-plane poles. Trials 3 and 4 represent 
second-order systems with real and complex poles respect- 
ively. Trial 5 represents a second-order system with 
complex poles and trial 6 is a third-order system with 

real poles of -10, -20 and -30. Only the five lowest-order 
coefficients of the estimates are given in Table 1 since 
the higher-order estimates become negligibly small. The 
results given are the averages of fifty consecutive esti- 


mates made for each system's coefficients. 


ae 


TABLE 1 


Brome r O Nrebrenoronenint 


Actual System coefficients = q, 


Estimated coefficients = d 


dq, wleugoruveda: d¢ = 0 
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VI. CONCLUSIONS 


The identification scheme presented in this investiga- 
tion does provide good estimates of the transfer function 
Sasa ions for an unknown system with no numerator 
dynamics. The estimated coefficients in Table l are seen 
to be within 50% of the actual coefficients for all cases 
except when the actual coefficient equals Zero. For these 
coefficients the estimates rapidly approach zero since the 
estimates are for higher order coefficients than those that 
exist in the actual system. The identification scheme 
thererore provrdes anmimearcarron of the order ™ol the “system. 

The need for a low pass filter on the input to the 
"unknown" system has not been substantiated. Trials 1 
through 6 of Table 1 indicate that the filter can be deleted 
with no disastrous effects. 

The feasibility of using the Fast Fourier Transform as 


a tool for system (plant) identification has been shown. 
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VIL. RECOMMENDATIONS 


This investigation has been unduly long and involved 
duewsto the requirement thal wall the modelime be done on a 
Gigital computer. Many of the parameters listed in Table 
Dl_of Appendix D required a large number of trial runs to 
gGetermine a best set of values for the overall simudlatginen.. 
These trials were necessary to aid in finding the minimum 
iteration rate and input signal sampling rate required to 
adequately simulate the continuous "unknown" plant since 
these rates are fundamental to the total simulation time 
(not to be confused with identification time). It is 
Suggested that any further investigation in this area be 
done on a Hybrid computer so that the continuous plant may 
be run onthe analog portion_of the computer and the iden= 
tification scheme on the digital portion. 

As a starting point the identification scheme requires 
that a guess be made as to the maximum number of coeffi- 
clients in the unknown plant's transfer function. <A major 
effect of this guess is that it determines the dimensions 
of matrix W of equation (15), Appendix C. Since the inverse 
Of Wis reduired it would appear that ats dimensions showils 
be kept as small as possible. However, it has been observeg 
that the dimensions of W that produce the most nearly cor- 
rect estimates of the transfer function coefficients are 
not the smallest (consistent with the number of coefficients 


known to exist) nor the largest (which depend on the amount 
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of computer memory reserved for W). The best dimensions 

for W are different even for different plants of the same 
order. It is suggested that any future investigation util- 
izing this identification scheme should involve a method 

of avoiding the inverse of W or attempt an analytic solution 


of wot. 


at 


APPENDiPS 4 
Thicwis cures Gia iee a Siicuh. oxen) 


The FFT is used to evaluate the complex Fourier trans- 
term of a sequence Of Complex numbers. It Y's a fimrcire 
discrete Fourier transform (DFT) with the restriction that 
the number of samples in the complex sequence is constrained 
to be an exact power of two. The power-of-two constraint 
on the number of samples allows implementation of an algor- 
1chm on 4a digital comouter for rapid computation. of sem 
Pinte WE. 

The mechanics of the FFT and its characteristics has 
been widely discussed in the literature [7, 8, 9, 10, ll, 
2. oe Ol Che Chatwaevctgmauices Will Ue presenuce 
IMIS - 

From Goldman [3], a function g(t) which exists only for 
-T a 


oi ae and which has a Fourier transform G(f) which exists 


only for -W <f<W can be expressed by the DFT pair 








teen 
TW J 
m ik nN 2'TW 
e(=-) == fF G(=) ¢€ (1) 
OW T -_ TW Ab 
Oo mn 
TW -j 
ike al m OTW 
G(r) = oy 2 &(sz7? = ee 
m=—-T'W 


where e(s,) = cGy, fermt = a 


—e sag 
G(r) = ewer wom cf 7 
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g(t) and G(f) are completely determined by their 2TW+1 
sample values EG) and Gi) respectively. 

The assumptions Jeading to (RP mancde( 2). thavee(t) exist 
iy Ole = <t <5 and G(f) exists only for -W<f<W, are 
physically unreadizable. For example j@eometder the case 


where g(t) is a section of a continuous function x(t) such 


mia c 
x(t) is Coatinueten tor ’—3@ <> <= (3) 
SG) = Gb) > yi 
y(t) = u(tt+s) - u(t-5) 


DUCE mse latem Ob ihn es eteon DE anonm nels 


Then Gp. = iis Cay <i —q ) da (4) 
2 sin 67 
and Y(t) = aie _aFT 5 (oy 


men) and therefore G(f) are continuous for =~ <f< =<, 


By observing (sampling) a section of x(t) such that 


eee) = x(t) for = <t <= the assumption leading to (1) 


and (2) are violated. 

ifva- y(t) ain (3 )aieehesen so that 2(v) = x(t)- yCC) 
has negligible value outside of the range a <t<5 and G(f) 
is negligible outside the range -W<f<W, then (1) and (2) 
@an be made to be approximately exact. y(t)'"s to satisfy 
the above include the Hanning and Hamming functions [4, 5] 


and the Parzen spectral window [14]. 


Zo 


For the FFT, (1) and (2) are normally expressed as 


N-1 4 emia 
g(m) = E G(nde s (6) 
n= 
o™mmn 
N= -j —— 
G(n) = = ze(m) & — (7) 
m= 
where g(m) = (sa) 
G(n) = G(m) 


Some authors put the weighting = on (6) rather than (7)8 
However, the FFT used in this investigation is from the 
TBM/360 Scientific Subroutine package (360A-CM-03X) and in 


this algorithm the = weighting is placed on (7). 
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ARREN DEX, B 
Discreve sOluumon Lor x(t) 


Given the vector differential equation (1) it is desired 


fowobtain a diisicrebeusoiution, fomeughs 


y= eg ereeecre sy Gre, is) 
T 
where x(t) = [x, (t) X(t) ---------- Xy(t ) J 
it 
GOS Ge) (Ge) =e Xy (6) 


Pesan OeN) onst aneemaTr i x 
Bis an (N,1) constant vector 


utc) mS a scalar 


The superscript T denotes the transpose of the vector. From 


Seana [2] theseontinuous time solution of (1) is 


A(t-t,) fe 
many =e O x(t,) +f ASP") Buz) at, (2) 
c 
@ 
ee) 
where e is the transition matrix o(t,t,) 


For a discrete solution of (2), where x(t) is sampled 


at intervals of time T, let 


Then 
x [ame = —° xen) + f Alt) T-tIe yay (3) 
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The usual procedure is to assume that v(t) = v(nT) for, 
nies + —(otlor. The imemememuation of (rT) am this: inves— 
tigation possesses characteristics which make it advisable 
CO GSalMplemunrs forcing funeumeonmmar as high a rate as 
possible. Since v(t) is assumed to be broadband noise, it 
is varying very rapidly with respect to time. Therefore T 
should be very small if v(t) is to be fairly constant over 
the sampling interval. Also it must be noted that in the 
simulation of the identification problem v(t) is represented 
by a sequence of numbers from a random number generator. 

The magnitude of the numbers in this sequence are indepen- 
denteof T in that if fi we members are produced whey ewer] 
have the same magnitudes whether they are produced at 
imperyals of T or T/e. Themimpli@ation then is that for 

any two consecutive numbers from the generator, such as v(i) 
and v(itl), those which are produced at T/2 second intervals 
represent a v(t) that has a higher rate of change with 
respect to time (and a higher frequency) than those that 

are produced atyT second intervals, i.e., 


v(itl)-v(i) 


T/d ; 


ier en 
a 


Therefore if the v's in (3) can be generated at intervals 
smaller than T, then the effective bandwidth of the noise 
can be made wider as the sampling interval of v(t) becomes 
smaller. Some averaging of the rapid fluctuations of v(t), 


weighted by the system parameters in the matrix A, will be 
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made by the integral in (3). The result is that v(t) can 
be sampled at a high rate (higher than 7) and the unforced 
response in (3) need only be computed every T seconds. To 


implement this approach let, 


(m+1 ) 


M~1 gi 
AG _ 
SiGe Tee aT") ae if ceA[(nt1)T-t]Bv(t)dt 
mee (Ben) T (4) 


A ie Beata Ns (stn) T <1 <(Sen)t 


Since B is @@constant vecter and vig) willbe agconstant 


Syer the interval of intec@earien , 


yer NT 
AT _ m 
x | Cail) Tt) = es ene Fe ii eA[(nt1)T-1 Jdt-Bev [Gtn)T] . 
m= U7 mi 
—— (5) 


To remove the parameter n from the integral in (5) let 


nT—=A 


= 
A = nT-T 
Then 
er 
AT ap Mol ON ar m 
x PCat) 'T ] =e 1 ge y -f € dA Bev [(qtn)T). 
m=O -mT 
M en) 


To remove the negative limits on the integral in (6) let 


_ mT 
XN = -T iM 
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Then (6) becomes 


Z T 
M-1 M -A(t+57 
sae eel cheer ee. (iO 


m=0 0 


AT AT 
E 


pe Gara Se a Bagnall ie 


Since the integral in (7) is not a function of m, then 


aL 


peinyy -AT -A 


M 
j@ <¢ She es Ta Bv[ (THN) T] (8) 
0 m=0 


AT 
E 


xpent Ley = x(n) Fe 


Wevevaluate the amcegral in @i on avdugiral comouveme 


-AT 2 . ° e ° ° e 
€ is expanded in an infinite series. 











te M 
= if ae CFL)! (9) 


ze 
m -AT 

EGuataen (9) 1s “tiemoelution ot {  « dt to be used. How- 
O 


ever, if A is known to be nonsingular, then 
al 
: -AT -1 —AT 


fi = dtm = A [ipe=a-! le 
0 


where I is the identity matrix of the same dimensions as A. 
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a-+ should exist except for those cases where the system 


characterized by (1) has a pole at S=0. However, to avoid 
compounding errors by taking an unnecessary matrix inverse 
on the computer, (9) is the solution of the integral to be 
used. 


From (8) ardm@(9)@and letting 4 = a 


AT AT e CS ee cara 
pach le) = 6 Ga ee \- a 8 -» Le ‘Bev(nT+md ) 
— — 4=0 (i4+1)! 5 
(ala) 
OR 
x[(nt1)T] = 6(T)x(nT) + P(T,M)v(nT,M) (12) 


where 


x(nT) is the (N,1) state vector of the plant at 
time t=nT 


$(T) is the (N,N) transition matrix Al 


T 1s the transition time between the system's states 
v(nT,M) is the (M,1) control vector | v(nT) 


v(nT+A ) 


v(nT+(M-1)A) 


/ 


r(T,M) is an (N,M) transfer matrix such that 


i er 5 
(-AX)™ | | 7° -AA* -2AA 7 (M=1) AX 


(Tr Z 


PCT .M) =a er 2 
i= 


ois 


APPEND A.C 
imnplenentavicmot=@vire’oystem Identificaticrselution 


A given system is represented by Figure Cl, where 


K 
Dic) =a! a. (jw)*. It is desired that the coefficients 
k=0 


of D(w) be identified through measurements on v(t) and x(t). 






wat.) 
V(w) 


Figure Ci. System Representation in the 
Frequency Domain 


Since 
dee X (w ) (1) 
D(w) Vlw) ? 
then 
X(w)D(w) = V(w) (2) 


Tf X(w) and V(w) are known for discrete values of w where 


Ww > We then 


X(w, )D(w, ) = V(w,) . Cy 


XCw, ) and VW, ) result trom apo lseaumen of the Fast 


Fourier Transform algorithm [Appendix A] to measurements 
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of S6CG)) axidwa.( bee eee re X(w, ) eine V(w, ) are assumed 
to be known. If D(w,) in (3) is replaced by D(w,)> the 


estimate of D(w.), then (3) becomes 
X(w,)D(w,) = V(w,) (4) 


Let r(w,) equal the residue, or error, due to the 


“~ 


estimate D(w,). 
r(w,) = X(w,)D(w,) - V(w,) (5) 


Then applying a least squares estimation procedure to 
- 
= ‘ ¥ 
Q : r(w,) r¥(w,) (6) 
* denotes the complex conjugate 


Q@ is to be minimized. 


From (5) and (6), 


i " x * 
Q = he [X(w,)D(w,)-VCw, )] [X(w, )D(w, )-VCw, )] 
(7) 
Let 

X(w, ) =a 

VCw, ) = V; 

a“ K oe 

D(w,) = me di) (ju, ) = D, 


oy 


iMieweesetving the partial of Q@ with respect to a coeffi- 


laa 


Cake ane On equal to zero, 

















6 I ar aD*, aD, . 
pee ey, CK Dea, 82 eee. (8) 
ee =O 9. wl od a © meee 
m ie m 
where 
aD, m 
ran (Jw, ) (9) 
ad 
m 
aD*, 5 “a Kk a m 
= — y dy (-Jw, ) = (-1) (jw, ) , (10) 
ad ad, k=0 


T 
Ny eaviitly? m i m 
z Vx ; 1) (jw, ) et , (Jo, ) (11) 


Since (j)" is common to both sides of (11) and is not a 


PUMC GaenmsOt 1. 


m 


L K mm _ 
EX xy EE (-1)™d, (jw, )* + (-1) "4, (jw, )“1 , 


T 
- _7)M. yx % mM 
z [(-1)". X#,V, + V8.X,] w,™ , (12) 
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where Dy has been replaced by its summation. The d's are 


not TUNCUTLONS, Ooi, cic 


Ca ees O, if m and k are not both even or odd 
= 2, if m and k are both even 


= -2 if m and k are both odd 


Therefore, for m = 0 to M and k = 0 to K the Yeft side of 


(12) can be expressed as 





QO ----- w(0,K) 
Gp) ee w(1,K) 
0) Qe = = 
TCs eee - el | le) 
nO w(M,K) 
where 
w(mk) = g* [(-2)™(-1)*] +w, OM) 
ae: 
= do 
oat 
qd5 
| d 
pes 


oe 


The right hand side of (12) can be expressed as 


| Re (X, * Ves ) 
~jim(X, eve ) 
Re ary 
2 (14) 
M 
My. #V.4V.4V. ¥X. : : 
toe 0 a Ge en, 7 





x = x 
where Re(X, V,) real part of X, V; 


He - e e x 
Im(X, V,) imaginary part of X, V, 
= Dk he — ello 

A ae a: 


T = the time span over which v(t) and x(t) are 


observed. 


Now to go from the foregoing general derivation to the 
particular case where the order of the system is N, let 
M= K=N. Then from (12), (13) and (14) and after simpli- 


fCac Oven  SOluiaiLOn 1 Om d becomes 
d=G Ww 2 (i 


where dq 
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By 
Eo 
q7+ = 
EN 
2) a iL 
ls 
Sy OT 
a me 
a pen a 
i 
0 
w(2,0) 
0 
ir + = as 
: w(4,0) 
w(N,0) 





aL H enleleignchise 
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w(N,1) 


-1l 





- | Z 
Z = ‘fala 
NG 
ie0 


Batu, a Jule for n even 


-Im(X, *V, )i° for n odd 


no =—0" to | 


The implementation of the identification scheme then 


takes on the form of Figure Ce. 


For example, if the order of the system N equals 2 and 


v(t) and x(t) are observed through a time window of T 


seconds then, 


dy nl 0 

Bt he Tv 
dy 0 Tea 
d., 0 0 -( 
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_a 
ii 
0 


ul 
x 
Re(X, V,) 
, . ¥ 
-Im(X, V, )4 
is Z| RetXx, *V ae 
; wD ails 
1=0 


v(t) ik x(t ) 





D(f) 
T 
-— He ee HE ANALOG 
Noa f 
are _ \ DIGITAL 
Va 
. 
PRODUCT CONJUGATE re, 
eat 
44 PRODUCT 
xX, 
SUMS AND 
PRODUCTS 
SUMS, 
PRODUCTS 
AND 
5 MATRIX 
en INVEPSE 
w-t 
ProDucT oF qvt wot » es 


Figure C2. Implementation of 
the Identification Scheme 


a 


Ne lec NO ID.« IB, 
Computer Program PLEST (PLANT ESTIMATION) 


The computer program is composed of the main program 
and subroutines GAUSS, DHARM, PLANT, RECIP and EA. GAUSS 
and DHARM are called from the IBM/360 Scientific Subroutine 
package (360A-CM-03X) and are not included in the program 
liste 

Subroutine GAUSS is used to obtain a sequence of real 
random numbers with a normal distribution and a selectable 
mean and standard deviation... The expected valine of Lhe 
power spectrum of the~output of GAUSS has been investigated 
by utilizing the Fast Fourier transform and has been found 
to approach a constant value over all frequencies. GAUSS 
is therefore used to approximate a source of "white" noise. 

Subroutine DHARM computes the Discrete Complex Fourier 
Transform of a sequence of numbers. DHARM is used to 
Obtain the Prequemcy=spectrum of the input and Burvput son 
the system (plant) for which the transfer function coef- 
fmerents are to De determined. 

Subroutine RECIP is used to obtain the inverse of a 
square matrix with dimensions of from (1,1) to(12,12). 


Subroutine EA computes either, 


5 AL en 
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Ole 


00 nN 
A 
2 ae 629) 
AEG i teeleaye 
where A isa sdtwere matrix 


When the sum of the magnitudes of all of the elements 


in the ne® 


ewe! the computation is stopped: 


term of (1) or (2) is less than or equal to 


Subroutine PLANT computes the 9% and [ matrices of 
equation (12) in Appendix B. The inputs to PLANT consist 
or eGivcvherwcnhe coefficiemes of a Cramsw#er Tuner on fora 
system or the A matrix and B vector from a vector differ- 


ential equation such as; 
x(t) = Ax(t)+Bv(t) (3) 


The MAIN program (PLEST) performs the normal bookkeeping 
functions required of any comprehensive program. It also 
utilizes subroutine PLANT so that a continuous system such 
as Figure Dl is converted to the discrete system in 
Figure D2. The MAIN program, subroutines DHARM and RECIP 
perform the operations required by Appendix C to identify 
whe system gnder test. 

The primary program parameters are identified in the 


table of Computer Program Data Cards (Table Dl). 
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Ge) 





P(t ,t,) -) v(t) 


a -& fool 
aS 
on 
Gr 
Gir 
(®) 
~UO 
ct iff 
1 
ct fe 
OK 


ame ef 2S SSS OOPS ele eee 


| 
+ Ct) 
| RCo tay = 
| —$tep | 
Ne ey ae | 

| : | 

[_ SvoltEM UNDER TEST = ae 
TO FOURIER TOV hOUn IER 
TRANSFORM TRANSFORM 

Ho oeee DIS Continuous Version of 


the System Under Test 
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NOISE GENERATOR 


u(T- n) (FROM SUBROUTINE 
7 GAUSS) 
> 8 _. 2 PSS SS] - 
4 1 woe CO) 
2 u(T) ‘ + W + DELAY 
| rn) CO a 
| A { | i. + | 
{ } 
| T.N ) ! o(T:N) | 
| ! ' | 
| 
3 r(T-(N~I)) | 
LN ucr-a-1) | 
LOW PASS FILTER = 
Ole(- ae le - nn MRR ES or 






% 


es 


r 
| 


la v(T-)7-(N-O 


| ToN-M-K 


OMe T 
HORAN ME che 
TRANSFORM 


Figure D2. Dis 


v(T.N) r(T-N) 


r(v-n-(M-1)) 


| | 
l | 
| | 
1 
| 





NOt: Sayer Cree Cr) T.NeM-K 


R=0 
aeletgu y(n, li) | 
from Equation (12) 
Appendix B 
Own fo 
POUR Tak 
TRANSFORM 


crete Version of the System under Test 


ta 


TABLE D1 


COMBUTERSEROGRAM DATAW CARDS 
Card Format Name Comments 


1 TY NJOBS Mie number of “sets Gf data cards 


to be operated on. 


2 10A8 Arbitrary comments to be read and 


written on the print-out. 


3 D10.4 VARN The variance of the random signal 
generated by the subroutine gauss 


(the mean has been preset to 0). 


4 DIG eh ae The time (T) between samples for 
the filver Amput . 


5 T4 N The number of successive inputs 
CoO the falter for OMe ourpum. 
The iteration time of the filter 
(TN) is determined by TN=T.N. 
Qk Ne< 12 


6 T4 M The number of successive inputs 
to the plant under test for one 
OUBDUT. new mveration tige oF 
the plant is determined by 
iN =e a. AO 


7 TY K For every Ke" iteration of the 


plant under test the plant input 
and output are sampled once for 
the anpuws to the FFT. “ies inpuce 
tO vCbemwur tT occur at inveryals=er 
ie ae Ne OS K Sake 
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Card Format Name Comments 
8 T4 NRUNS The number of times the plant 
coefficients will be deter- 
mined. The resulting coeffi- 
CLEMES are wl homoaverarte wel 
the results from NRUNS esti- 


mates of the coefficients. 


9 T4 LTRMS A lower-limit guess on the 
number of termnssin the plane 


Transfer funce von. 
1 < LTRMS <S)aeeis 


alge’ T4 NTRMS An upper-limit guess on the 
number of terms san eine ola 
tener teat Wierenu alice iar. 
1 See e 


‘Lal 14 CSURSUL The number of iterations the 
filter-plant combination will 
go through to allow it to 


redchma steady sstavemc oma lr erm. 


eZ T4 M(1) oM(1) samples of the plant 
IMDUt Vande GulpUrearerr oO. le 
used on each run (each NRUNS 
from card 8) to estimate the 


transfer function coefficients. 
SG) a0 


ahs} 2 tee The method by which the fil- 
Perewiik be eh aigaei Cr zeieee T 
ee one Cransier f£unemien 
COeiiacients are CoO game entered. 
if [Fer =c tae Desvom row or 
the A matrix and each element 
Sir une Bb veCcror is Co be 


entered. 
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Card “Formac Name Comments 


14 


as 


NODE: 


i). 1G 


(1) 16GB 


(1) 16- 


Cay 


(lf 


of a 


NOTES 


ee) 164A 


ihe ND Pie ieee The exponent of highest power 
of S in the denominator of 
the Ss Miver teansier Tuner won 
0 < NDR esi 


iy VerreT The exponent of highest». power 
of S inthe numerator of the 
filter Cranster fume Garon. 
OR NEA ss DRE LT 


if PLT = eee tou(2). 168 


blo. per The coefficient of the highest 


power ef S in the denominagvom 


DLO Dees) 


Dioe mo peCNDPEETe) “The goerrictent of S~ amneene 


transfer function denominator. 


DLC SUN GL } The coefficient of the highest 
power of S in the transfer 


runmCE LOM MuMetau er. 


Did, eeaee2) 


D10.4 UN(NPFLT+1) ine mC Oc iiivine ken ie cOg 39 in the 


transfer function numeravor. 


1G ee, fexee ccc) ite) 


ple, SCN Dr fy) ae The element A(NDPFLT,1) in 
Cie wine heer A mat rink. 
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Card Format Name 


(2 sie 


=) 2c. 
2) 174 


ce) 178 


ce) 17- 


18 


D1O0.4 A(NDPFLT,2) 


AOS BUN ONDIE ESE SIN IDE a Lg) 


DIO. 4+ SEG) 


Die. S22) 


D1O0.4 B(NDPFLT) 


Ih Jeol 


Comments 


The element A(NDPFLT,NDPFLT) 
inet. SPilter AY mate x. 


The elements 1) Of time 


Filter Bevecpor 


The elemen@ B(NDPFLYT) of 
the filter B vector: 


The method by which the 
plant wili@pe charac. — 
1Zed. Soeame scomme mir cms 
for Camas 


Cards 19 through 22 (forgthe plant) ame similar ao 
cards 14 through 17 (for the filter). 


PORE : 


Caio. 

5 

13 

14 

AWS 
el aro k 
Cr) 7 a 


to omit the filter enter the following data 


Data 

0001 

ul 

00 

00 

TY. COgoD+00 
Ihe QIOLOKO IRE O10 


Repeat Data Cards 2 through 22 NJOBS (from Card 1) times. 
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